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Abstract 

The first part of the report summarizes research carried out under 
this grant and now in the published literature; included are basic quantum 
detection and estimation theory, applications to optics, photon counting, 
and filtering theory. The second part describes recent work on the 
restoration of degraded optical images received at photoelectrically emissive 
surfaces; the data used by the method are the numbers of electrons ejected 
from various parts of the surface. 
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I. Survey of Published Results 


The research accomplished under this grant has in large part been 

published in the papers listed on pages 9-11. This final report will be 

mostly devoted to describing our recent work on the restoration of degraded 

optical images received at photoelectrically emissive surfaces, but first 

we recapitulate our previous investigations. In particular we call 

* 

attention to our review paper in Progress in Optics (§1, 3) , which 
summarizes the elements of quantum detection theory and its application 
to the detection of light sources and the estimation of their parameters. 
Earlier reviews (§1, 1,2) described the state of the theory at the time, 
particularly with regard to its application to communications. 

The function of the receiver in a communication system is to decide 
which of a set of known signals has been transmitted during each element 
of a uniform sequence of intervals; it issues its decisions as a sequence 
of symbols representing its best version of the transmitted message. In 
a radar the purpose of the receiver is to decide whether the echo from a 
target is present in its input and, if so, to estimate the range and 
velocity of the target. In optical communication and radar systems the 
information on which such decisions and estimates are based is embodied 
in the field at the aperture of the receiver, where it is corrupted by 
random background radiation. Because that field is governed by the laws 
of quantum mechanics, the methods of classical decision and estimation 
theories cannot be applied, and those theories must be reformulated to 


•k 

This notation refers to paper 3 of §1 of the bibliography on pages 9-11, 
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take into account quantum-mechanical restrictions on measurability 
of the fields (§1, 3). 

In quantum detection, the decisions made by the receiver are 

considered as choices among a set of density operators p., p , p , 

12 M 

each describing the aperture field in the presence of one of the possible 
transmitted signals. The density operators act in the Hilbert space <30?^ 
of the field. The choices are based on measurements of the field and 


calculations involving their outcomes , and the entire procedure can be 
viewed as applying to the field a resolution of the identity operator 1 


in into a set of positive-definite Hermitian operators, 11^, n^, 



i 


n i + II 2 + +II M = l 


such that the probability that the receiver decides that the j-th signal 
was sent when really the k-th was sent is the trace Tr(p R II^). The optimum 
receiver will apply that resolution of the identity for which the average 
probability of error is minimum. It has been shown that the operators 
11^, n 2 , ..., n M attaining minimum error probability need not commute, but 
that when they do not, an ancillary apparatus can in principle be coupled 
with the receiver in such a way that the same minimum error probability 
can be attained by measuring commuting projection operators in the Hilbert 
space of the combination (§1, 5). 

When the receiver estimates parameters of a received optical signal, 
such as its arrival time and Doppler shift, errors are introduced by the 
background radiation and by the stochastic quantum-mechanical nature of 
the signal itself. Lower bounds on the mean-square errors of such estimates 
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can be set by quantum-mechanical counterparts of the Cramer-Rao inequality 
of classical statistics. Two forms of the inequality are known, and 
sometimes the one, sometimes the other yields the superior bounds (§1, 5). 

Quantum detection theory was originally formulated in terms of an 
ideal receiver consisting of a lossless box or cavity into which the 
incident light was admitted by opening an aperture during an observation 
interval (0,T). Measurements were then to be made on the electromagnetic 
field inside the cavity at a later time t > T in order to decide whether 
it contained a component attributable to the source or object sought, or 
in order to estimate certain parameters, such as the radiance or direction 
of an optical source. Ordinarily an optical instrument bases its decisions 
and estimates about a source on the electromagnetic field at its aperture, 
for instance by focusing the incident light onto a photographic plate or 
a counter whose response is presented to an observer. We wanted, therefore, 
to eliminate that artifice of the lossless cavity and to express detection 
and estimation strategies directly in terms of the aperture field. This 
was done first for threshold detectors (§11, 1) and estimators (ill, 2). 
Optimum detectors and estimators , as well as Cramir-Rao bounds on mean-square 
errors of estimates of parameters, must in quantum mechanics be based on 
the density operators of the aperture field, and by decomposing that field 
into spatio-temporal modes it has been possible to write down the necessary 
density operators in a broad class of problems involving coherent or 
naturally incoherent light (ill, 3; §1, 3, §4 and §5). 

The resolution of close point sources of incoherent light could be 
conveniently studied by means of this modal decomposition of the aperture 



field. Two formulations were treated by quantum detection theory. In 
the first, an observer is to decide whether two point sources of known 
radiant power and position are present, or whether only a single source 
having twice the radiant power and located at the midpoint is present. 

In the second formulation either one or the other of two point sources 
may be radiating during the observation interval, as in a binary 
communication system in which 0's and l's are sent by turning on one 
source or the other. In both formulations the density operators do not 
commute. The basic eigenvalue equation of quantum detection theory was 
written as an integral equation by expressing the density operators in the 
coherent-state representation, and by assuming the absence of background 
radiation it was possible to solve the integral equations for both 
formulations in closed form and calculate the error probabilities incurred by 
the optimum decision strategies (§1, 4). 

The information of interest in a scene emitting or reflecting natural 
light is embodied in its radiance distribution as a function of position 
and frequency, and a camera or a telescope is essentially an instrument 
for estimating that radiance distribution, which the laws of light 
propagation translate into the spatio-temporal coherence function of the 
field at the aperture of the instrument. The restoration of images degraded 
by diffraction, aberrations, and atmospheric turbulence can be viewed as 
a problem of estimating a comprehensive set of parameters of the aperture 
field or of the radiant source producing it. Treatment of this problem by 
conventional linear estimation methods, assuming additive Gaussian noise, 
is often unjustified. Within the framework of classical optics it was 


shown how the random fluctuations of the light fields cause large errors 
in estimates of the radiance of an object plane at points closer together 
than half the Rayleigh resolution distance (§11, 4). Conditions were 
derived under which a formula of the Shannon type holds for the information 
transfer from an incoherently radiating object plane (§11, 5). An analysis 
was made of a system that estimates the radiance distribution of an object 
plane by focusing it on a dense array of photoelectric detectors and 
processing by linear means the numbers of photoelectrons emitted by each. 
Equations of the Wiener type were obtained for the optimum linear processor, 
with the noise term embodying the classical fluctuations of the object and 
background light fields, and the quantum fluctuations of the photoemissive 
process (§11,6). 

The ideal quantum receiver of a coherent signal of random phase counts 
the number of photons in a single mode of the aperture field, or in a mode 
of the lossless cavity suitably matched to the signal. The number of 
photons has a Laguerre distribution in the presence of a signal and a 
Bose distribution in its absence. Detection and error probabilities for 
such an optimum receiver have been calculated (§111, 1). In practice 
the signal will be passed through a narrowband filter in order to cut out 
as much background light as possible and will then be focused on a photon 
counter; detection and error probabilities for such a system have also 
been calculated (§111, 2). 

For many detectors, and particularly for the likelihood-ratio detectors 
pre'scribed by detection theory, it is not possible to calculate analytically 
the probability distributions of the statistics on which are based decisions 



about the presence or absence of a signal, but one can work out their 
moment-generating functions , which are the Laplace transforms of the 
distributions . Methods for determining the cumulative distribution of a 
statistic from its moment-generating function are therefore of interest, 
and two have been studied, one involving an expansion in terms of Laguerre 
functions (§IV, 1), the other the evaluation of the inverse Laplace 
transform by the method of steepest descent (§IV, 3). 

These techniques must be applied in analyzing the performance of a 

•k 

photoelectric detector proposed some time ago. An observer is to decide 
whether a certain image is present or not amid background light incident 
on a photoelectrically emissive surface divided like a mosaic into many 
small, isolated elements. The decisions are based on the numbers of 
electrons emitted from each element, and these have Poisson distributions. 
The likelihood ratio is formed from the numbers and compared with a 
decision level chosen, for instance, to yield a pre-assigned false-alarm 
probability. The moment-generating functions of the logarithm of the 
likelihood ratio can be expressed in terms of the illuminances caused by 
the background and by the image sought, and by finding their inverse 
Laplace transforms by the method of steepest descent, false-alarm and 
detection probabilities have been calculated. This detector has been 
compared with the threshold detector^ and with one that simply counts the 
total number of photoelectrons emitted from a circle concentric with the 


C. W. Helstrom, "Detection and resolution of optical signals". 
Trans . IEEE , vol. IT-10, pp. 275-287 (October, 1964). 
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expected image and having optimum radius; an image with a Gaussian 
distribution of illuminance was postulated. It was found that the optimum 
likelihood-ratio processing of the photoelectrons yields only slightly 
greater detection probability than that of the simple counter (§111, 3). 

In a thesis completed under this grant (§V, 5), Dr. C. L. Rino, now 
of Stanford Research Institute, developed and analyzed methods for solving 
the convolutional integral equation describing, inter alia , the linear 
filtering undergone by an image in passing through an isoplanatic optical 
system. It is this integral equation that must .be solved in restoring 
degraded optical images by linear estimation. Most attention was given 
to one-dimensional problems. Use of the discrete Fourier transform was 
investigated (§V, 3,4), and the absence of data beyond the endpoints of 
the data interval was taken care of by using causal Wiener filters derived 
for the semi-infinite interval. The extrapolation of tbe solution beyond 
the ends of the data interval was studied by means of expansions in terms 
of prolate-spheroidal wavefunctions , and it was shown how the properties 
of those functions explain the large errors inherent in such extrapolations 
(§V, 1). These functions were also used to treat the detection and 
estimation of bandlimited signals (§V, 2). 
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II. RESTORATION OF DEGRADED OPTICAL IMAGES* 

1. Maximum-Likelihood Image Restoration 

The restoration of degraded optical images has in the past been treated 
largely as linear filtering of a spatial random process in the presence of additive 
Gaussian noise. ^ ^It is only under special conditions that this is an accurate 
model of the image. Although the electromagnetic fields of the light from the ob- 
ject and the background combine additively, it is not the field of the object light 
that is of interest, but the radiance distribution of the object or scene, which is 
related to the covariance function of the field. Furthermore, what is measured at 
the image plane of an optical system is the illuminance of the incident light, which 
is also a quadratic functional of the field. When photographic film is used, the 
noise depends on the illuminances due to object and background in a complicated way. 
The simple additivity of the Gaussian model is seldom realistic. 

In an endeavor to study image restoration in terms of a more realistic, yet 
mathematically tractable model, we have postulated that the incident light is fo- 
cused on a photoelectrically emissive surface, which is divided like a mosaic into 
a number N of small, insulated regions A^ of area AA. During an observation in- 
terval of duration T the numbers n^ of photoelectrons ejected from each of these 

elements A. are recorded and constitute the data on the basis of which an estimate 
x 

of the undegraded image is formed. 6 The system is illustrated in Fig. 1.1. 

Let the illuminance that would appear at the point x of the image plane, 
were the undistorted geometrical image present, be B(x); this is proportional to 

* 

The research described in this part was carried out in collaboration with 
Y.-M. Hong and is described in his doctoral thesis, "Optical Signal Processing: 
Poisson Image Restoration and Shearing Interferometry", Dep't. of Applied Physics 
& Information Science, University of California, San Diego; September, 1973. 
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Radiance B 



Fig. 1.1 Image restoration system. The elements s of area AA emit numbers 
n^ of electrons that are processed to yield estimates B • B.. + b 
of samples of the true image of the radiance distribution B. 0 » 
object plane, A = aperture, I = image plane. A narrowband spectral 
filter for object and background light is not shown. 
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the radiance distribution of the object plane with its dimensions scaled ac- 
cording to the magnification factor of the optical system. The actual illumin- 
ance at point x is 


I(x) = / K(x,u) B(u)d u. 




( 1 . 1 ) 


where K(x,u) is the incoherent point spread function of the optical system and 
the transmission medium combined. It accounts for turbulent distortions of 
the refractive index, diffraction at the aperture of the observing instrument, 
aberrations in lenses, and relative motion of object and instrument. We suppose 
the kernel K(x,u) to be known. 

Let _x^ be the coordinates of the center of the i-th element of the image 
plane, which element we suppose small enough that I(x) is effectively constant 
over it. Then the mean number of phofoelectrons emitted by the i-th element A 


n i = “ I ^ i >» 


i 1,2,.. ,n. 


( 1 . 2 ) 


with 


a=ThkJ 


n(v)$(v) (hv) dv, 


( 1 . 3 ) 


where n(v) is the quantum efficiency of the surface at frequency v, $(v) is the 
relative spectral density of the light, 


$(v)dv = 1, 


and hv is the energy of a quantum of light at frequency v. When as we assume 
here the duration T is much greater than the reciprocal W ^ of the bandwidth of 
the light, WT >> 1, the numbers n_^ of emitted photoelectrons have Poisson dis- 
tributions and are statistically independent, ^ 
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rl n n i 

Pr J"i "h| ' [J, ( -"i ) 


(1.4) 


We suppose that 


B(u) = B q + b(u). 


(1.5) 


where B^ is known and b(u) is to be estimated at a number M of uniformly 
spaced points u., in the intervals between which B(u) changes in a smooth 

* ~i >x> 

manner. With 


= f 


u) B n d u 
o ~ 


( 1 . 6 ) 


assumed independent of x, we put for the mean values of the data n^. 


n i = a(I o + n i )s 


(1.7) 


and 


where 


M 

”i -£ Vj 

j=l 


i — 1,2, ... ,N, 


( 1 . 8 ) 


K. = K(x. ,u.) AA. 


(1.9) 


Thus in this approximation 
M 

= 


E K. . B_. 
ij 0 

j=l 


( 1 . 10 ) 


The data are the N numbers n^ distributed as in (1.4); the estimanda . are the 
sample values b^. of the deviation of the geometrical image from a mean level 

V 

In order to apply the principle of maximum likelihood, some statistical 
description must be formulated .for the class of expected true images B(u). We 
consider them as spatial Gaussian random processes with mean zero and auto- 
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covariance function 

' £ [ B( »i >B( ii2>] - B 0 2 " i [ b fai) b( J!2 ) ) o.ii) 

The images are spatially stationary, and the width of <p(u) represents the size 
of typical details in the object plane; the ratio 
V(0) 

C (1.12) 

B q + <?(0) 

specifies a mean-square contrast. The Fourier transform of the autocovariance 
function (p( u) describes the distribution of spatial frequencies in objects of 
the class. Although the actual objects will not in general resemble what we 
think a Gaussian stochastic process looks like., a method that works well for 
an ensemble of such processes should usually be effective in restoring an object 
of similar structure. The structure of a scene can seldom be specified in terms 
less rudimentary than the average contrast and the size of typical details. The 
form chosen for the autocovariance function in (1.11) represents the demands we 
are placing on our restoration scheme. If (p (u) is taken to be very narrow, we 
are forcing the method to try to reproduce fine details, and we can expect the 
resulting image to be noisy. If <f> (u) is too broad, the resulting image will be 
in error because details have been smeared out. 
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where 


q (N) ■v xp (■ * gk -ik b j \) 

is the prior p.d.f. of the estimanda, Pr ^jn^j |^j | j 
distributions given in (1.4), and 


(1.14) 

is the product of Poisson 


Pr (M) '/ •■•/ q (N) Pr (hl|i b il) db i-- db M 


is the unconditional distribution of the data. In (1.14) the y., are the elements 

jk 

of a matrix y = </> ^ that is the inverse of the covariance matrix cp whose elements 


are 


nj = ^ ( iii ' ~j ^ ’ 


(1.15) 


Cq in (1.14) is a normalizing constant. Maximizing the conditional p.d.f. in 
(1.13) is equivalent to minimizing its negative logarithm 


N 


M M 


f (S i " n ± ln V + ? EE VA- 


(1.16) 


i=l 


j=l k=l 


I h; 


from which terms independent of the estimanda |bj| have been dropped. The 
constraints on the minimization are embodied in eqs. (1.7) and (1.8), which take 
the distortion of the image into account. 

If the numbers n^ are large and the contrast is small, the distributions 
of the nAs are approximately Gaussian with mean values Tu given by (1.7) and 
variances n^ ~ a Iq. The quantity to be minimized is then, instead of (1.16), 

f* f|b,n - i (ay- 1 £ (n i - y 2 

M M 


+ 7 EE VA- 

j=l k=l 


(1.17) 
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which must also be taken with (1.7) and (1.8). The minimizing values of the 


estimanda are the solution of the linear equations 


N 

1 

i=l 


V ( I 0 P jk + K ik K ij) b k (n i " “V K ij 

k=l 


N 

I 

i=l 


(1.18) 


they provide the usual least-squares estimate, which can be written in matrix 
form as 


A / T \-l T, 

b = I n n + aK K K (n - 
~ y 0~ j ~ ~ 


*V> 


■ ( T Ai * a Z^) (a - aI oi) 


(1.19) 


where n is an N-dimensional column vector of the data, 1 a column vector of 

N l's, and 1^ is the MxM identity matrix. The matrix K = | l^jl I is NxM and 
T 

K is its transpose. This least-squares estimate forms the starting point of 
our search for the maximum-likelihood estimate, which minimizes F^jbjjj in (1.16) 
Having neither prior knowledge nor preconceptions about the class of 
images to be restored is equivalent statistically to attributing an infinite 
variance to the estimanda jb^j. The second term of (1.16) then disappears, 
and one minimizes only 

N 


F "(! b ji) (“i - "i ,e -i) 

i=l 


( 1 . 20 ) 


in combination with the constraints in (1.7) and (1.8). For the least-squares 
estimate this leads to the simple "inverse filtering" specified by the equations 


b = (a K T K ) “ 1 K T (n - al l) 


(i.?d 



We shall find that the image restorations based on prior assumptions about the 
class of true images, through (1.16), are seldom much more accurate than those 
obtained, by maximizing (1.20), without utilizing such assumptions. 
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2. The Method of Steepest Descents 


The most successful method we tried for computing the maximum-likelihood 
estimate ^ of the true image was the method of steepest descents applied to the 
minimization of the function F ^ j b ^ | ^ in eq. (1.16).® Let a given set of trial 

A 0 

values of the estimates be denoted by the superscript . We expand the function 

F ^ | b ^ | ^ in an M-dimensional Taylor series about the point jb_.^j, 

'(M- ft)" (v4 


k=l 


+ 4 


M M 

w 


2 Zj/j \ 3b, 3b 

ii i 'km 

k=l m=l 




(2.1) 


and treat this as a quadratic form to be minimized. The next trial point lies 
in the direction of the gradient from Jb^, 


b. = b. - A g., 
J J J 


( 2 . 2 ) 


where 


hi 

i=l i k=l 


y jk b k°’ j < 2 * 3 > 


are the components of the gradient vector at point ]3 , in which we must put, 
from (1.7) and (1.8), 


”1 ■ “ ( x o + h 0 )’ 

h ‘It Vj°’ 


j-i 


(2.4) 


(2.5) 
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The parameter X determines how far along the gradient one proceeds, and it is 
determined by substituting (2.2) into (2.1) and minimizing with respect to X, 


M / H M 

a *> /ee 

/ k=l m=l 


( 2 . 6 ) 


_ / 3 F \ 

\m \ 9b, 9b / 

' k m ' 


^km + 


n . K ,K 
l ik im 

( V *! 0 ) 2 


The quadratic form in the, denominator of X can be written 


(2.7) 


where 


M M M . M M n, r 

EE ^m^k^m EE %av + E 


k=l m=l 


.-t 


k=l m=l 


i? (W? 


C ik s k* 


( 2 . 8 ) 


(2.9) 


When the trial point changes from Jj to b by (2.2), in (2.5) changes from 
n i ° to ”^ r i * T ^ ie onl y matrix that needs to be stored in the computer 
is y = <P , the inverse covariance matrix, and for the covariance matrix^ we 
used, y was very simple. The rest of the operations in calculating X and the 
new values of |b_.| and | n j | can be carried out with one-dimensional arrays. 

Because the true image must be positive, values of b^ that sank below 
-Bq at any stage of the iteration were set equal to -B^-t-e, where e was a small 
positive constant. The iterations were stopped when the ratio 


t 'v; 


fell below 10 -3 , where |b^°'| were the values at the beginning, |b^ those 


at the end of an iteration. 
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3. Generation of Trial Objects and Images 

The method of steepest descents was tried out on the restoration of 
one-dimensional objects distorted by motion blur and diffraction. Two types 
of objects were used, a deterministic one composed of Gaussian peaks and random 
ones generated by an autoregressive process. The former had the form 

b^ = 90 | exp f”-(i-5) 2 /4^j + exp |"-( 


i-9> 2 /4 


J + exp J^— (i— 16) 2 / 4 J | . (3.1) 


Its peaks are two units in width and located at coordinates 5, 9, and 16 in a 
total range of twenty points . 

The random object was generated by the formula 

b/_^2 — r i 2, ...» M— 1, (3.2) 


starting from an initial value b^ that was a Gaussian random variate with mean 

2 

zero and variance o^ . The z/s were independent Gaussian random variates with 

2 2 - 

mean zero and variance a. (1-r ) , as a result of which all the image variates 

b 

2 

. were Gaussian with mean zero and variance a, , but with a correlation matrix 
' b 

given by 


2 i-j 

% = "b r 


(3.3) 


The inverse V of this matrix is especially simple, 
diagonal elements are 

y ll = y MM = °b /( ' 1 ~ V ) ’ 

V J ii = a b 2 (l+r 2 )/(l-r 2 ) , i ^ 1 or M, 

^i,i + l = y i + l,i = - °b’ 2 r/(1 ’ r2) - 


Its diagonal and first off- 


M 
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All other elements of the matrix y vanish. 

The Gaussian variates z^ were generated from pairs of uniformly dis- 
tributed random variates and R 2 in the interval (0,1) by means of the 
tr ans f o rmation 

z = a z (-2 In R 1 ) 1 ^ 2 cos (2 ttR 2 ), = (Var z) 1 ^ 2 (3.5) 

Values of b^ generated by (3.2) that fell below -Bq were set equal to -B^ in 
order to produce an object with nonnegative radiance B^ + b . A typical object 
generated by this process is illustrated by the heavy line in Fig. 3.1. Its 
parameters are B^ = 3, r = 0.6, = 25. 

When the image is degraded by relative motion between object and re- 
cording instrument, the point spread function K(x,u) will be constant for a 
range of values of x about ij and zero elsewhere. Sampling then produces an 
NxM matrix K of the form 

rw , 

K ±j = D _1 , i = j , j+1, . . • , j+D-1, j = 1,2 M, (3.6) 


where N = M+D-l, and D is an integer proportional to the relative velocity of 
object and instrument; we took D=4. Although (3.6) appears to violate (1.10), 
we can assume that the actual object extends to the left and right of the range 
of sampled values and has in those margins a constant radiance B^; the remaining 
analysis, with Iq = B^, is unaltered. 

For diffraction the point-spread matrix was taken as 


K.. = C j^sin 2 W(i-j)J / (i-j) 2 = e.^ | i- j | <J , 

= 0 , 1 i-j I > J, 


(3.7) 
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25 


Again we suppose that the object has wide enough margins with b(u) = 0, in 
which its radiance is constant and equal to Bq = Iq. Here W is related to the 
width of the aperture, which is taken to be a slit perpendicular to the one 
dimension of the object. 

The data n., representing the numbers of photoelectrons ejected by 

X / 

the incident light from the several elements of the image surface, were generated 
as Poisson-distributed random variates with mean values 


n i " a + V* 


i. = > K. .b., i = 1, . . . ,N. 
L-t xj j 


(3.8) 


Let 1 = n. be a typical such mean value. A pseudorandom-number program produces 
a uniformly distributed random variate y in (0,1), which is multiplied by e\ 

I 

The sums \ 



j-0 

are accumulated, starting with s_ = 1, until s just exceeds ye\ whereupon 

u n 

the accumulation stops and the datum n^ is set equal to the final integer n. 
The probability that n^, takes on a value n is then, with s ^ = 0, 

Pr |n t = nj = Pr js^^ ye X < s n J = X n e ~ X /nl, 


as required. 

For mean values ru exceeding 25 we approximated the Poisson distribution 

as a Gaussian with mean value n^ and variance n^; a Gaussian random variate was 

— 1/2 

generated as in (3.5) with o z = n^ , n^ was added, and the result was rounded 



to the nearest integer. 

The object shown as the heavy line in Fig. 3.1 was blurred by applying 

the matrix K in (3.6) with D=4, and Poisson variates n. were generated from the 
~ i 

resulting mean values n^ obtained as in (1.8), with a = 10. The n^'s were 
then divided by a to reduce them to the same scale as the object and plotted 
on Fig. 3.1 as the thin wavy line, which represents the recorded data on which 


the restoration was based. 
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4. The Computer Simulations 

The image data generated as we have just described were subjected 
to three image-restoration procedures and the results were compared. As a 
simple measure of the quality of a restoration we chose the relative average 
squared error, in percent, defined by 

1 00 xf: (Vj) 2 /e b j 2 • <4 - i) 

j=i j=i 

where b. are the known illuminance values of the true image and b. are the 
2 2 

estimates. 

The first restoration procedure was simple inverse linear filtering, 
which corresponds to the least-squares estimate in the limit <p (0) -»•“>; no 
prior knowledge of the class of objects is presumed. It is determined from 
eq. (1.21). The second procedure was the method of linear least-squares es- 
timation, as given by eq. (1.19). The resulting estimate^ served as the 
starting point for seeking the nonlinear maximum- likelihood estimate, and for 
this the method of steepest descents described in Section 2 was applied. For 
all three procedures any estimates b^ lying below -Bq were set equal to -Bq 
before the average squared error & was calculated. 

Figure 4.1 shows the results of applying these three procedures to 

restoring the image depicted in Fig. 3.1. As the signal- to-noise ratio 
2 

aa^ /Bq is large in this case, the linear least-squares estimate is not much 
different from the one obtained by inverse filtering as in (1.21); the 
relative percentage errors are 10.1% and 12.5%, respectively. The maximum- 
likelihood estimate, on the other hand, yields $ = 2.53% and lies much closer 
to the true image than either of the others. 




When establishing an image-restoration system it is often difficult to 

judge what covariance function <p(u) most aptly describes the class of true 

images to be restored. The robustness of the system to the use of inappropriate 

covariances is therefore of interest. In order to investigate it, we have tried 

2 : 

our methods with values of the object variance a, and the correlation coefficient 

b ; 

r in (3.3) and (3.4) different from those with which the object in Fig. 3.1 was 
generated. The resulting error percentages & are plotted in Fig. 4.2. Inverse 
filtering is independent of the object covariance. The maximum-likelihood method 
is relatively insensitive to the value of the correlation coefficient r until it 
exceeds about 0.9, after which the method breaks down and yields larger errors 
than the linear least-squares estimator. It may be that a value of r too close 

to 1 imposes too much correlation on the estimated radiance values and prevents 

! ' 

the search procedure from accommodating to the data. 

The absence of any prior assumptions about the class of images being 


restored is equivalent to the assignment of infinite variance to the prior p.d.f. 


M 


(!-,!) 


to 


in (1.14) and eliminates the second term from the function F 
be minimized, eq. (1.16). The steepest-descent procedure minimizing (1.20) is 
therefore simpler and cheaper. Fig. 4.3 shows the resulting restoration of the 
image in Fig* 3.1. Its relative squared error amounts to 2.69%, only slightly 


larger than for the restoration that used the known values of c^ ! and r and is 


depicted in Fig. 4.1. 

The variability of our restorations was studied by starting with the 
same distorted image, that is, with the same mean values n., and generating 
independent sets of Poisson-distributed data n.^ by starting the random number 
procedure with different "initializing constants." The first example used an 
autoregressive object with 0 b =15, r = 0.2, I Q = 3, and a = 10, distorted by 
motion blur introduced by the matrix K of (3.6) with D=4. Table I lists the 



SQUARED ERROR (PERCENTAGE) 


Fig. 4. 


I 



O.OIO.l 0.3 0.5 0.6 0.7 0.9 0.99 


r 

✓ 

Percentage squared error in restoring Fig. 3.1 as function of 
assumed correlation coefficient r; curves are indexed by assumed 
standard deviation o^. Solid lines-linear estimate; dashed lines 
maximum-likelihood estimate. 
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Table I 



Q 

O' 

II 

15, r = 0.2, I Q = 3, 

a = 10 




Relative Squared Errors 

(percent) 

* 


°b " 15 

, r = 0.2 

a b •** 



linear 

nonlinear 

linear 

nonlinear 

1. 

5.53 

1.80 

8.81 

0.862 

2. 

2.13 

2.46 

1.02 

1.08 

3. 

4.90 

2.93 

4.76 

2.21 

4. 

5.95 

4.28 

8.20 

2.44 

5. 

7.61 

5.20 

5.17 

3.43 

* 

6. 

5.83 

5.79 

3.43 

3.58 

7. 

8.49 

4.23 

8.25 

3.68 

8. 

6.02 

5.64 

3.97 

3.69 

9. ' 

5.49 

3.67 

9.77 

4.76 

10. 

9.84 

6.22 

8.69 

5.17 

11. 

7.84 

6.36 

7.24 

5.63 

12. 

12.59 

9.39 

12.36 

9.56 

mean 

6.85 

4.83 

6.81 

3.84 


standard 2.67 2.07 3.18 

deviation 

■k 

Illustrated in Figs. 4.4 and 4.5 
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resulting squared errors; the left-hand columns are for a restoration procedure 

using the correct prior variance = 15 and the correct correlation coefficient 

r=0.2. The right-hand columns are for i.e., no prior knowledge. The 

linear and the maximum-likelihood estimators were used; for a,-*® the linear 

b 

estimator corresponds to inverse filtering as in (1.21). The starred line 
refers to the restorations depicted in Figs. 4.4 and 4.5. 

In ten out of twelve of these trials the smaller relative error & was 
v obtained for 0^-*°°, i.e., by discarding all prior knowledge of the class of ob- 
jects. In two cases the linear estimator was more effective than the maximum- 
likelihood estimator. In general, however, the maximum-likelihood estimator 
was rather better than the linear estimator, and as the figures for standard 
deviation at the bottom of the table demonstrate, the former is the more con- 
sistent. 

A second set of trials with a different object of the same type yielded 
the squared errors listed in Table II; the starred case is illustrated in Figs. 
4.6 and 4.7. Here there was a greater variability, and in all but one case the 
maximum-likelihood estimate was closer than the linear one to the true image. 
Omitting the prior information led to closer restorations in seven out of the 
twelve trials. Again the maximum-likelihood estimator was the more consistent, 
and on the average considerably better than the linear estimator. 

The Gaussian "deterministic" object of eq. (3.1) and its motion-blurred 
image are shown in Fig. 4.8; again the matrix K in (3.6) was used with D=4. 

The dip between the two left-hand peaks has been wiped out. Data were 
generated from this object by taking = 1^ = 5 and a = 3. It was then 






Table II 


a b = 15, r = 0.3, I Q = 5, a = 10 


Relative Squared Errors (percent) 



V 15 > 

r = 0.3 

a b'*“ 



linear 

> nonlinear 

linear 

nonlinear 

1 . 

7.91 

6.70 

7.74 

5.42 

2. 

9.69 

7.07 , 

8.11 

5.70 

3. 

14.93 

8.84 

11.82 

6.59 

4. 

13.57 

7.23 

15.79 

6.90 

5. 

9.16 

5.10 

11.55 

7.51 

6. 

13.94 

8.12 

18.44 

9.73 

* 7. 

19.62 

11.75 

19.91 

10.17 

8. 

9.51 

6.86 

9.94 

10.74 

9. 

29.62 

8.62 

47.06 

10.99 

10. 

16.04 

14.55 

16.81 

14.27 

11. 

23.18 

18.07 

24.70 

17.02 . 

12. 

22.19 

12.78 

25.73 

17.05 

mean 

15.78 

9.64 

18.13 

10.17 

standard 

deviation 

6.70 

3.85 

10.90 

4.11 


Illustrated in Figs. 4.6 and 4.7 











restored under the assumption that it belongs to the class of objects gen- 
erated by the autoregressive process of (3.2) with a, = 45 and r = 0.5. The 

b 

results of the maximum- likelihood and linear least-squaresestimates are shown 

by the crosses and circles in Fig. 4.9. In Figs. 4.10 and 4.11 we plot the 

relative squared errors obtained when different values of a, and r were used 

b 

in the restoration procedure. For most of those values the maximum-likelihood 

estimate is. the more accurate. It is also less sensitive to the choice of the 

values of o^. and r than the linear least-squares estimate. 

Diffraction was the second mode of image distortion considered. Figure 

4.12 shows an autoregressive object generated with the values a, = 15 and r = 0.2 

b 

and distorted by the diffraction kernel of eq. (3.7) with W=1.2. The small central 

peak was eliminated by the diffraction, and we were unable to perceive it in any 

of the restorations. Data were generated from this image with 1^ = 3 and a = 10. 

The restoration procedure assumed the values a = 20 and r = 0.2, producing the 

points shown in Fig. 4.13. The three large peaks stand out in both restorations; 

the other two did not survive. In Fig. 4.14 are shown the percentage squared 

errors attained for various assumed values of the standard deviation a. and the 

b 

correlation coefficient r. The method of inverse filtering failed completely for 

this object. Discarding prior knowledge (o, = °°) , the maximum-likelihood estimate 

b 

yielded a relative squared error of 24.52%, which is roughly equal to the smallest 

value a ttai ned by any of th e linear methods. ' 

In Fig. 4.15, finally, we show the relative squared errors attained in 

restoring a diffracted Gaussian object of the form in eq. (3.1). The value of W 

in (3.7) was taken equal to 0.8, and I n = 10, a = 10. In the restoration a, was 

U b 

taken equal to 15, and various values of r were tried. Over most of the 
range the maximum-likelihood estimator is much superior. 
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Fig. 4.10 Relative squared errors in restoring Fig. 4.8 for various assumed 
values of standard deviation r = 0.5. Solid line-linear least 
squares; dashed line — maximum likelihood. 
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Fig. 4.11 Relative squared errors in restoring Fig. 4.8 for various assumed 
values of correlation coefficient r ;a^ = 45. Solid line-linear 
least squares; dashed line-maximum likelihood. 





Inverse filtering 48 80 % 


d = 00 


24.52 % 



Relative squared error (percent) as function of correlation 
coefficient r for restoring image in Fig. 4-12. Curves are 
indexed by the assumed standard deviation o fe . Solid lines- 
linear least squares; dashed lines-maximum likelihood. 
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Fig. 4.15 Relative squared error (percent) in restoring image of eq. (3.1) 
after diffraction with W » 0.8 as function of assumed correlation 
coefficient r. I Q » 10, a - 10, a b » 15. 


48 


5. Conjugate Gradient Method 


When the method of steepest descents nears the minimum value of the 


function F 


(hi)- 


it has a tendency to oscillate between two points in the 


space of the estimanda, and this behavior protracts convergence. In an ef- 
fort to find a faster computing routine, we turned to the method of conjugate 


gradients. We applied it only to the function F 


"(hi) 


in (1.20), which is the 


negative logarithm of the likelihood function in the absence of prior knowledge 
about the class of objects. This method was unsuccessful. 

Again a quadratic approximation to the function to be minimized is used, 
as in (2.1). The new trial point is derived from the original one at any stage 
by shifting in the direction d, which is one of the so-called" conjugate directions". 


b = b + Ad. 


(5.1) 


The value of A is again chosen to minimize the new value of F", 


A = - d-g / d T Ad, 


(5.2) 


where g is the gradient vector whose components are given in (2.3), and A is the 
matrix whose elements are a^ m given in (2.7); in both cases we have set the in- 
verse covariance matrix p equal to 0. The first direction d is taken as -g; 
thereafter the direction d for a stage is given in terms of the direction d' for 
the previous stage by 


d = yd' - g, 


(5.3) 


where 


T' / , T , 

y = gg/g g 


(5.4) 
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g' being the gradient vector at the previous stage. After M such stages, d 
~ ~ 

is once again set equal to -g, the displacement of steepest descent, and the 
cycle is repeated. 

The difficulties with this method apparently arose from the positivity 
constraint, which required our setting b^. = -B^ + e whenever b^. fell below -B^. 
Equation (5.1) then required us to change the corresponding component of the 


vector d' to 


(- b o v)A 


* A before calculating the new directional from (5.3). 

This occasionally caused the next move to be in an unfavorable direction, and 
the method did not converge. The images on which we tried it had high contrast, 
so that there were several adjacent points with zero illuminance, B^ + b^ = 0. 
For an object of low contrast, with b^ never falling below -B^, the conjugate- 
gradient method doubtless would have been superior, but for such objects the 
least-squares solution in (1.19) is generally adequate. The nonlinearity of the 
positivity constraint and the logarithmic singularity of the function F"||b^||to 
be minimized are too far from compatibility with the basic assumption of the 
conjugate-gradient method — that the function is a quadratic form — for it to be 
successful with objects of high contrast. 


6. Conclusion 


In the great majority of cases tried, the maximum-likelihood estimate 
of the true image was superior to that calculated by the linear least-squares 
method. The minimization of the negative logarithmic likelihood function by 
the method of steepest descents appears to be the most efficient procedure. 
The prior information about the class of images being restored had in most of 
our trials little influence on the outcome, and minimization of the function 



( 6 . 1 ) 


( 6 . 2 ) 


b . 
J 



(6.3) 


provided good estimates of the true images in most trials. This insensitivity 
to prior assumptions is gratifying, for the proper statistical description of a 
class of images is difficult to formulate. 

This maximum-likelihood estimate is most appropriate when the image has 
a high contrast, with few or no photoelectrons being ejected from substantial 
areas of the picture. It is under such circumstances that the usual linear model 
of image formation with additive Gaussian noise is least accurate. The function 
in (6.1) takes the Poisson distributions of the data n^ into account, and the 
mechanism by which the image was distorted is embodied in the constraint (6.2). 
The third constraint (6.3) forces the restored image to have everywhere a non- 
negative illuminance. There remains to be investigated how the computation can 
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be made more efficient when applied to pictures furnishing large numbers of 
data. 
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Figure Captions 


Fig. 1.1 

Fig. 3.1 
Fig. 4.1 
Fig. 4.2 

Fig. 4.3 

Fig. 4.4 

Fig. 4.5 

Fig. 4.6 

Fig. 4.7 

Fig. 4.8 
Fig. 4.9 

Fig. 4.10 
Fig. 4.11 


Image restoration system. The elements s of area AA emit numbers 
n_^ of electrons that are processed to yield estimates B = Bq + b 
of samples of the true image of the radiance distribution B. 0 = 
object plane, A = aperture, I = image plane. A narrowband spectral 
filter for object and background light is not shown. 

True image (thick line) and the noisy degraded image (thin line) 
after blurring by relative motion. Iq = 3, a = 10, a = 25, r = 0.6 

True image of Fig. 3.1 and its restorations: o = linear estimate, 

A = inverse-filtering estimate, O = maximum- likelihood estimate. 

Percentage squared error in restoring Fig. 3.1 as function of 
assumed correlation coefficient r; curves are indexed by assumed 
standard deviation cr^. Solid lines-linear estimate; dashed lines — 
maximum-likelihood estimate. 

Restoration of Fig. 3.1 without prior assumptions (o^ = °°) . o = 
linear inverse filtering, x = maximum-likelihood. 

True image (thick line) and data (thin line) after distortion by 
relative motion, a = 15, r= 0.2, Iq = 3, ct = 10, D= 4. 

Restoration of image in Fig. 4.4: x = maximum-likelihood; o = 
linear inverse filtering, = ■». 

True image (thick line) and data (thin line) after distortion by 
relative motion, o = 15, r= 0.3, Iq = 5, a= 10, D = 4. 

Restoration of image in Fig. 4.6: x = maximum likelihood; o = 

linear inverse filtering, c. = °°. 

b 

Deterministic object of eq. (3.1) and its motion-blurred image; D = 

Restoration of image in Fig. 4.8: Iq = 5, a = 3, = 45, r = 0.5. 

d = maximum- likelihood estimate, o = linear least-squares estimate. 

t 

Relative squared errors in restoring Fig. 478 for various assumed 
values of standard deviation a^; r = 0.5. Solid line-linear least 
squares; dashed line — maximum likelihood. 

Relative squared errors in restoring Fig. 4.8 for various assumed 
values of correlation coefficient r ;a^ = 45. Solid line-linear 
least squares; dashed line-maximum likelihood. 

Diffraction of an autoregressive image: a = 15, r = 0.2, I = 3, 

i rt i o ^ 


Fig. 4.12 



Fig. 4.13 Restoration of image in Fig. 4.12, taking a, = 20, r = 0.2, 

W = 1.2. o - linear least squares; 4- maximum likelihood. 

Fig. 4.14 Relative squared error (percent) as function of correlation 
coefficient r for restoring image in Fig. 4.12. Curves are 
indexed by the assumed standard deviation . Solid lines — 
linear least squares; dashed lines-maximum likelihood. 

Fig. 4.15 Relative squared error (percent) in restoring image of eq. (3.1) 
after diffraction with W = 0.8 as function of assumed correlation 
coefficient r. I- = 10, a = 10, o, = 15. 



